#include "rand.h"
#include <cmath>
#include <cstdlib>

namespace randy
{
    // The following code was taken from: https://gist.github.com/1354700
    // Original author: DRMacIver
    // Merged with some code from:
    // And from http://clerc.maurice.free.fr/Maths/index.htm
    double random_normal();


    void random_sphere(int n, double *data) 
    {
        double tot = 0;
        double x;

        for(int i = 0; i < n; i++) {
            x = random_normal();
            //x = mtrand();
            tot += x * x;

            data[i] = x;
        }
  
        double norm = std::sqrt(tot);
        double rd = 1.0 / static_cast<double>(n);
        double t = random_uniform();
        double radius = std::pow(t, rd);
        for(int i = 0; i < n; i++){
            data[i] = radius * data[i] / norm;
        }
    }

    double random_uniform()
    {
        return ((double)std::rand()) / ((double)RAND_MAX);
    }

    double random_normal()
    {
        double u, v, s;

        // Box-Muller transformation (polar coordinates)
        do {
            u = 2 * random_uniform() - 1;
            v = 2 * random_uniform() - 1;
            //u = 2 * mtrand() - 1;
            //v = 2 * mtrand() - 1;
            s = u * u + v * v;
        } while(s >= 1);

        return u * std::sqrt(-2 * std::log(s) / s);
    }
    
}
